% 参数设置
m = 1000;
n = 500;
k = 300;
num_trials = 10;  % 重复试验次数

% 预分配结果存储
times_qril = zeros(num_trials, 1);
times_pinv = zeros(num_trials, 1);
errors_qril = zeros(num_trials, 1);
errors_pinv = zeros(num_trials, 1);
diff_pinv = zeros(num_trials, 1);

fprintf('===== 伪逆算法性能测试 =====\n');
fprintf('矩阵大小: %d×%d, 系数维度: %d×%d\n', m, n, n, k);
fprintf('试验次数: %d\n', num_trials);
fprintf('--------------------------------------------\n');
fprintf('试验\tQRIL时间(ms)\tpinv时间(ms)\tQRIL误差\tpinv误差\t与pinv差异\n');
fprintf('--------------------------------------------\n');

for trial = 1:num_trials
    % 生成随机数据
    A = rand(m, n);
    beta_std = rand(n, k);
    T_true = A * beta_std;
    
    % 计算最小维度
    d = min(m, n);
    T = eye(d);
    
    % 测试QRIL算法
    tic;
    if m > n
        Beta_qril = QRIL_inverse.QR_schmidt(A', T);
        Beta_qril = Beta_qril';
    else
        Beta_qril = QRIL_inverse.QR_schmidt(A, T);
    end
    times_qril(trial) = toc * 1000;  % 转换为毫秒
    
    % 测试MATLAB pinv
    tic;
    Beta_pinv = pinv(A);
    times_pinv(trial) = toc * 1000;  % 转换为毫秒
    
    % 计算误差
    beta_recon_qril = Beta_qril * T_true;
    beta_recon_pinv = Beta_pinv * T_true;
    
    errors_qril(trial) = norm(beta_recon_qril - beta_std, "fro");
    errors_pinv(trial) = norm(beta_recon_pinv - beta_std, "fro");
    diff_pinv(trial) = norm(Beta_qril - Beta_pinv, "fro") / norm(Beta_pinv, "fro");
    
    % 输出当前试验结果
    fprintf('%d\t%.4f\t\t%.4f\t\t%.2e\t%.2e\t%.2e\n', ...
            trial, times_qril(trial), times_pinv(trial), ...
            errors_qril(trial), errors_pinv(trial), diff_pinv(trial));
end

% 计算平均结果
avg_time_qril = mean(times_qril);
avg_time_pinv = mean(times_pinv);
avg_error_qril = mean(errors_qril);
avg_error_pinv = mean(errors_pinv);
avg_diff_pinv = mean(diff_pinv);

% 输出汇总结果
fprintf('--------------------------------------------\n');
fprintf('平均\t%.4f\t\t%.4f\t\t%.2e\t%.2e\t%.2e\n', ...
        avg_time_qril, avg_time_pinv, avg_error_qril, avg_error_pinv, avg_diff_pinv);

% 可视化结果
figure;

% 时间对比
subplot(2, 2, 1);
bar([times_qril, times_pinv]);
title('计算时间对比');
xlabel('试验编号');
ylabel('时间 (ms)');
legend('QRIL', 'pinv', 'Location', 'northwest');
grid on;

% 误差对比
subplot(2, 2, 2);
semilogy(1:num_trials, errors_qril, 'b-o', 'LineWidth', 1.5);
hold on;
semilogy(1:num_trials, errors_pinv, 'r-s', 'LineWidth', 1.5);
title('系数重建误差');
xlabel('试验编号');
ylabel('Frobenius范数误差');
legend('QRIL', 'pinv', 'Location', 'northwest');
grid on;

% 与pinv差异
subplot(2, 2, 3);
semilogy(1:num_trials, diff_pinv, 'm-d', 'LineWidth', 1.5);
title('与pinv结果的相对差异');
xlabel('试验编号');
ylabel('相对Frobenius范数');
grid on;

% 时间分布
subplot(2, 2, 4);
boxplot([times_qril, times_pinv], 'Labels', {'QRIL', 'pinv'});
title('计算时间分布');
ylabel('时间 (ms)');
grid on;

% 伪逆矩阵可视化对比
figure;
subplot(1, 3, 1);
imagesc(Beta_pinv);
title('MATLAB pinv结果');
colorbar;

subplot(1, 3, 2);
imagesc(Beta_qril);
title('QRIL算法结果');
colorbar;

subplot(1, 3, 3);
imagesc(abs(Beta_qril - Beta_pinv));
title('差异绝对值');
colorbar;

% 重建系数对比
figure;
subplot(1, 3, 1);
imagesc(beta_std);
title('原始系数');
colorbar;

subplot(1, 3, 2);
imagesc(beta_recon_qril);
title('QRIL重建系数');
colorbar;

subplot(1, 3, 3);
imagesc(beta_recon_pinv);
title('pinv重建系数');
colorbar;